Background

Having protein alternatives is important for sustainable scaling of agriculture¹ to meet the growing food demand. Given the amount of coastal area in many countries, understanding ways to use those resources is important. Marine aquaculture has the potential to play a key role in increasing food production, boosting economic growth in coastal areas, and help maintain clean waterways².

Marine aquaculture potential is dictated by many variables including ship traffic, dissolved oxygen, bottom depth³. Understanding what coastal areas have potential for different species can guide the implementation of marine aquaculture.

The Exclusive Economic Zones (EEZ) represent areas from coastlines that countries have resource jurisdiction over, and could target for developing aquaculture. In this project, I am determing the locations in West Coast EEZ that are suitable for a variety of oyster aquaculture. The areas need to fit these specific conditions:

sea surface temperature: 11-30°C
depth: 0-70 meters below sea level

Goals

Determine the suitability area in the Exclusive Economic Zones of the West Coast
Illustrate what EEZ have the greatest percentage and area for suitable marine aquaculture
Create workflow to determine suitable marine aquaculture area for other species

Datasets

The stars and terra packages will be useful for handling the raster data.

1. Sea Surface Temperature

This analysis will characterize sea surface temperature within the regions with sea surface temperature (SST) from 2008 to 2012. The data is originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.

This data is stored in the sst_tifs folder: - average_annual_sst_2008.tif
- average_annual_sst_2009.tif
- average_annual_sst_2010.tif
- average_annual_sst_2011.tif
- average_annual_sst_2012.tif

2. Bathymetry

To characterize the depth of the ocean we will use the General Bathymetric Chart of the Oceans (GEBCO).⁴

read in bathymetry raster depth.tif for depth data

3. Exclusive Economic Zones

Data on the Exclusive Economic Zones off the west coast of the US are from Marineregions.org.

load in the wc_regions_clean.shp file to access the EEZ shapefiles

Load libraries

# load libraries
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(here)

## here() starts at /Users/lunacatalan/Documents/dev/eds223/final_repos/oyster-aquaculture-suitability

library(stars) # for .tif files

## Loading required package: abind

library(terra)

## terra 1.7.55
## 
## Attaching package: 'terra'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(leaflet)
library(patchwork)

## 
## Attaching package: 'patchwork'
## 
## The following object is masked from 'package:terra':
## 
##     area

library(shiny)
library(htmltools)

Read in the Data

Read in the EEZ and depth data:

# load in data for west coast regions
wc_eez <- st_read(here("data/wc_regions_clean.shp")) %>% 
  rename(ID = rgn_id)

## Reading layer `wc_regions_clean' from data source 
##   `/Users/lunacatalan/Documents/dev/eds223/final_repos/oyster-aquaculture-suitability/data/wc_regions_clean.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
## Geodetic CRS:  WGS 84

# read in bathymetry ocean depth data
depth <- read_stars(here("data/depth.tif")) %>% 
  rast()

Since there are multiple .tifs for the sea surface temperature data, I made a list and made a raster stack instead of individually reading in the files and then stacking them. Simplified!

# read in annual average sea surface temperatures
file_names <- list.files("data/sst_tifs/",
                         full.names = TRUE)

# stack the rasters together
sst_stack <- rast(file_names)
summary(sst_stack)

##  average_annual_sst_2008 average_annual_sst_2009 average_annual_sst_2010
##  Min.   :280.4           Min.   :279.4           Min.   :280.2          
##  1st Qu.:285.3           1st Qu.:285.7           1st Qu.:285.7          
##  Median :287.1           Median :287.7           Median :287.4          
##  Mean   :287.1           Mean   :287.7           Mean   :287.5          
##  3rd Qu.:289.1           3rd Qu.:289.7           3rd Qu.:289.4          
##  Max.   :301.4           Max.   :301.5           Max.   :301.0          
##  NA's   :42125           NA's   :42204           NA's   :42229          
##  average_annual_sst_2011 average_annual_sst_2012
##  Min.   :278.9           Min.   :278.5          
##  1st Qu.:285.5           1st Qu.:285.6          
##  Median :287.0           Median :287.0          
##  Mean   :287.1           Mean   :287.2          
##  3rd Qu.:288.8           3rd Qu.:289.0          
##  Max.   :307.3           Max.   :310.2          
##  NA's   :42071           NA's   :42034

class(sst_stack)

## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"

Check the crs for the data

Check the crs of the data before continuing so that we dont run into trouble when we try to make them interact:

# do the crs match?
st_crs(wc_eez) == st_crs(sst_stack) # no they do not
st_crs(wc_eez) == st_crs(depth) # yes

The crs for the West Coast data and the Sea Surface temp rasters don’t match. We have to re-project so that we can work on them together.

Reproject

To make sure that we can plot the sea surface temperature with the rest of the data we reproject to the crs of the west coast eez raster.

# reproject sst_stack to same crs as wc
sst_stack <- project(sst_stack, 
                     crs(wc_eez))

# check if the new object has new crs
summary(sst_stack)

##  average_annual_sst_2008 average_annual_sst_2009 average_annual_sst_2010
##  Min.   :280.4           Min.   :279.4           Min.   :280.2          
##  1st Qu.:285.3           1st Qu.:285.7           1st Qu.:285.7          
##  Median :287.1           Median :287.7           Median :287.4          
##  Mean   :287.1           Mean   :287.7           Mean   :287.5          
##  3rd Qu.:289.1           3rd Qu.:289.7           3rd Qu.:289.4          
##  Max.   :301.4           Max.   :301.4           Max.   :300.9          
##  NA's   :42125           NA's   :42206           NA's   :42229          
##  average_annual_sst_2011 average_annual_sst_2012
##  Min.   :278.9           Min.   :278.5          
##  1st Qu.:285.5           1st Qu.:285.6          
##  Median :287.0           Median :287.0          
##  Mean   :287.1           Mean   :287.2          
##  3rd Qu.:288.8           3rd Qu.:289.0          
##  Max.   :307.3           Max.   :309.9          
##  NA's   :42071           NA's   :42034

Prepare the depth and sea surface temp data

To use them together we need to make sure that the resolutions and extents of the objects are the same.

Make raster with just the mean sea surface temp over the 2008-2012 time period.

The temperature data is in Kelvin. To avoid having to transform the inputs of temperature parameters into Kelvin, we can change the means to Celsius bu substracting by 273.15 degrees.

# mean from stack
sst_mean_K <- terra::mean(sst_stack, 
                          na.rm = TRUE)

# convert values into C by subtracting by 273.15 
sst_mean <- sst_mean_K - 273.15

sst_mean

## class       : SpatRaster 
## dimensions  : 480, 408, 1  (nrow, ncol, nlyr)
## resolution  : 0.04165905, 0.04165905  (x, y)
## extent      : -131.9848, -114.9879, 29.99208, 49.98842  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :      mean 
## min value   :  5.041986 
## max value   : 32.720087

Get depth data with the same resolution as the sea surface temp mean and crop to the same extent as the sea surface temp raster.

Use resample and crop from the terra package to match the resolutions and extents of the sst and depth data:

# resample depth to match resolution of sst
depth_resample <- resample(depth, sst_mean,
                           method = "near")

# crop the depth raster to the sst raster
depth_crop <- crop(depth_resample, sst_mean)

Check the resolutions and extents of the depth and sst rasters:

# compare what the extents are visually 
plot(depth)
plot(sst_mean)
plot(depth_crop)

# check if the crs match
st_crs(depth_crop) == st_crs(sst_mean)

# look at the resolutions
depth_crop
sst_mean

Determine which locations are suitable for Oyster aquaculture using depth and temperature parameters.

Reclassify temperature raster:

I am using the sea surface temperature that oysters can suitably grow in: 11-30 deg C - recalssify so that all values outside of the 11-30 degree range are NA

# sea surface temperature: 11-30&deg C
rcl_temp <- matrix(c(11, 30, 1, # from EQUAL TO 11 to EQUAL TO 30, set the value to 1
                     -Inf, 11, NA, # from infinity to 11, set the value to 0
                     30, Inf, NA), # from greater than 30 to infinity, set the value to 0
                   ncol = 3, byrow = TRUE)

# reclassify
rcl_sst <- classify(sst_mean, 
                    rcl = rcl_temp)

Reclassify depth raster:

I am using the depth range that oysters can suitably grow of 0-70 meters below sea level. This translates to -70 and 0 meters. Any values outside of this range will be NA.

# depth: 0-70 meters below sea level
rcl_depth <- matrix(c(-70, 0, 1, # from EQUAL TO -70 to EQUAL TO 0, set the value to 1
                      -Inf, -70, NA, # from infinity to less than -70, set the value to 0
                      0, Inf, NA), # from greater than 0 to infinity, set the value to 0
                    ncol = 3, byrow = TRUE)

# reclassify
rcl_depth <- classify(depth_crop, 
                      rcl = rcl_depth)

Combine the reclassified depth and temp layers to find areas that are suitable for oysters:

Since the values in each raster are not the values of 1 (if suitable) and NA (if not suitable), we can perform map algebra and multiply the layers. Anywhere that there are NA’s will return an NA value, and only the pixels that have a value of 1 in both layers will return a 1. Therefore, the resulting raster will only have 1’s in the areas that match the suitability parameters in the depth AND the temp conditions.

# create multiplying function : layers that have 0 will multiply to 0 - only the places where it equals 1 in both will be 1 after multiplying 
rast_mult <- function(layer1, layer2) {
  (layer1*layer2)
}

#create stack of the reclassified data
depth_sst_stack <- c(rcl_sst, rcl_depth)

# apply the function to each layer in the stack
depth_sst <- lapp(depth_sst_stack, # stacked data
                  fun = rast_mult) # function

#rename layer as suitable
names(depth_sst) <- "suitable"

# check if it got renamed
depth_sst

## class       : SpatRaster 
## dimensions  : 480, 408, 1  (nrow, ncol, nlyr)
## resolution  : 0.04165905, 0.04165905  (x, y)
## extent      : -131.9848, -114.9879, 29.99208, 49.98842  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : suitable 
## min value   :        1 
## max value   :        1

Determine suitability in the EEZ to investigate which areas have the most potential for oyster marine aquaculture.

To work with the EEZ data we can rasterize it using the raster we created above with suitable areas.

This allows us to mask the EEZ data with the depth_sst data to identify what areas within the EEZ are suitable. To find the area that these pixels represent we can use the expanse function that sums the raster cells that are 1’s and not NA’s.

# rasterize the eez areas
id_rast = rasterize(vect(wc_eez), # create a SpatVector to make into a raster
                    depth_sst, 
                    field = "rgn") # select layer in database

# cells that are in the right spot
eez_suitable <- mask(id_rast, depth_sst)

# Compute the area covered by polygons or for all raster cells that are not NA in km
eez_zonal <- expanse(eez_suitable, 
                     unit = "km", # output is km2 because its an area
                     byValue = TRUE)

We can then calculate the percent of each EEZ that is suitable for oyster aquaculture, and comapre it to the total area in each EEZ.

# find the percentage of each zone that is suitable
eez <- wc_eez %>% 
  mutate(area = eez_zonal$area,
         percent = (area/area_km2)*100) # find percent

Suitability area visualization!

This gives a guide for regions to understand the productivity/economic potential of the area, and quantifying the total area is good for understanding priority for investing in marine oyster aquaculture.

Map for the Percent Suitable Area in the Economic zones:

Create basemap:

# basemap 
basemap <- leaflet() %>% 
  addProviderTiles(
    "Esri.WorldImagery",
    group = "Esri.WorldImagery"
  ) %>%
  # add a layers control
  addLayersControl(
    baseGroups = c("Esri.WorldImagery")
  )

Leaflet map with basemap of percent suitable areas:

pal <- colorNumeric(
  c("#E1F5C4", "#B3E0A6FF", "#7DC370FF", "#59A253FF","#368747FF"),
  # colors depend on the count variable
  domain = eez$percent,
)


map_percent <- basemap %>% 
  addPolygons(data = eez,
              # set the color of the polygon
              color = ~pal(percent),
              # set the fill opacity
              fillOpacity = 0.9
  ) %>%
  # add a legend
  addLegend(
    data = eez,
    pal = pal,
    values = ~percent,
    position = "bottomleft",
    title = "Percent Suitable \nArea:",
    opacity = 0.9
  ) %>% 
  addControl("Percent Suitable Area for Oysters In West Coast EEZ",
             position = "bottomright",
  )

map_percent

Map for the Total Suitable Area in the Economic zones:

pal <- colorNumeric(
  c("#F4D166FF", "#F8AF50FF", "#F38C30FF", "#EB6C1CFF","#CB4D22FF"),
  # colors depend on the count variable
  domain = eez$area,
)


map_total <- basemap %>% 
  addPolygons(data = eez,
              # set the color of the polygon
              color = ~pal(area),
              # set the fill opacity
              fillOpacity = 0.9
  ) %>%
  # add a legend
  addLegend(
    data = eez,
    pal = pal,
    values = ~area,
    position = "bottomleft",
    title = "Total Suitable \nArea:",
    opacity = 0.9
  ) %>% 
  addControl("Total Suitable Area for Oysters In West Coast EEZ",
             position = "bottomright",
  )

map_total

Develop a function to automate the processing of determining suitable area for marine aquaculture of other species on the west coast.

Re-set the west coast eez zones file:

wc_eez <- st_read(here("data/wc_regions_clean.shp")) %>% 
  rename(ID = rgn_id)

## Reading layer `wc_regions_clean' from data source 
##   `/Users/lunacatalan/Documents/dev/eds223/final_repos/oyster-aquaculture-suitability/data/wc_regions_clean.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
## Geodetic CRS:  WGS 84

The function `species_maps` :

Takes in 5 parameters: species, temp_low, temp_high, depth_low, and depth_high. - Will produce two visualizations of maps that show percent rnak and total area rank for EEZ on the west coast of the US.

species_maps <- function(species, temp_low, temp_high, depth_low, depth_high) {
  
  # sea surface temperature:deg C
  rcl_temp <- matrix(c(temp_low, temp_high, 1,
                       -Inf, temp_low, NA,
                       temp_high, Inf, NA),
                     ncol = 3, byrow = TRUE)
  
  # reclassify
  rcl_sst <- classify(sst_mean, 
                      rcl = rcl_temp)
  
  # depth meters below sea level
  rcl_depth <- matrix(c(depth_low, depth_high, 1, 
                      -Inf, depth_low, NA,
                      depth_high, Inf, NA),
                      ncol = 3, byrow = TRUE)
  
  # reclassify
  rcl_depth <- classify(depth_crop,
                        rcl = rcl_depth)
  
  rast_mult <- function(layer1, layer2) {return(layer1*layer2)}
  
  #create stack of the reclassified data
  depth_sst_stack <- c(rcl_sst, rcl_depth)
  
  # apply the function to each layer in the stack
  depth_sst <- lapp(depth_sst_stack, # stacked data
                    fun = rast_mult) 
  
  # rasterize the eez areas
  id_rast = rasterize(vect(wc_eez), # create a SpatVector to make into a raster
                    depth_sst, 
                    field = "rgn") # select layer in database

  # cells that are in the right spot
  eez_suitable <- mask(id_rast, depth_sst)

  # Compute the area covered by polygons or for all raster cells that are not NA in km
  eez_zonal <- expanse(eez_suitable, 
                     unit = "km", # output is km2 because its an area
                     byValue = TRUE) %>% 
    rename(rgn = value)

  
  # find the percentage of each zone that is suitable
  eez <- wc_eez %>% 
    right_join(eez_zonal, by = "rgn", copy = TRUE) %>% 
    mutate(percent_area = (area/area_km2)*100) # find percent
  
  pal <- colorNumeric(
  c("#E1F5C4", "#B3E0A6FF", "#7DC370FF", "#59A253FF","#368747FF"),
  # colors depend on the count variable
  domain = eez$percent,
  )
  
  percent <- leaflet(eez) %>%
    addPolygons(data = eez,
    # set the color of the polygon
    color = ~pal(percent_area),
    # set the fill opacity
    fillOpacity = 0.9) %>%
    addLegend(
      data = eez,
      pal = pal,
      values = ~area,
      position = "bottomleft",
      title = "Total Suitable \nArea:",
      opacity = 0.9) %>% 
    addControl(paste("% Suitable Area <br> for", species), # updates with species input,
             position = "bottomright",
             ) %>% 
    addScaleBar() %>% # add scalebar
    addTiles() # add basemap
  
  pal <- colorNumeric(
  c("#F4D166FF", "#F8AF50FF", "#F38C30FF", "#EB6C1CFF","#CB4D22FF"),
  # colors depend on the count variable
  domain = eez$area,
  )

map_total <- basemap %>% 
  addPolygons(data = eez,
    # set the color of the polygon
    color = ~pal(area),
    # set the fill opacity
    fillOpacity = 0.9
  ) %>%
  # add a legend
  addLegend(
    data = eez,
    pal = pal,
    values = ~area,
    position = "bottomleft",
    title = "Total Suitable \nArea:",
    opacity = 0.9
  ) %>% 
    addControl(paste("Total Suitable Area <br> for", species), # title that updates with species input
             position = "bottomright") %>% 
  addScaleBar() %>% 
  addTiles()

# put the leaflet maps next to each other
  map_grid <- 
    tagList(
      tags$table(width = "100%",
        tags$tr(
          tags$td(percent),
          tags$td(map_total)
        )
        )
      )
  
  # show the maps
  return(map_grid)
  
}

Testing the function on a new range of temps and depths!

# testing the function 
species_maps("Clam", 6, 12, -200, -30)

Fisheries, N. (2022, September 15). Aquaculture supports a sustainable earth. https://www.fisheries.noaa.gov/feature-story/aquaculture-supports-sustainable-earth ↩︎
Aragão, Cláudia, et al. “Alternative Proteins for Fish Diets: Implications beyond Growth.” Animals : An Open Access Journal from MDPI, U.S. National Library of Medicine, 7 May 2022, www.ncbi.nlm.nih.gov/pmc/articles/PMC9103129/.↩︎
GEBCO Compilation Group (2022) GEBCO_2022 Grid (doi:10.5285/e0f0bb80-ab44-2739-e053-6c86abc0289c).↩︎
GEBCO Compilation Group (2022) GEBCO_2022 Grid (doi:10.5285/e0f0bb80-ab44-2739-e053-6c86abc0289c).↩︎

Oyster Aquaculture Suitability in EEZ on the West Coast

Luna Herschenfeld-Catalan

2022-12-3

Background

Goals

Datasets

1. Sea Surface Temperature

2. Bathymetry

3. Exclusive Economic Zones

Load libraries

Read in the Data

Check the crs for the data

Reproject

Prepare the depth and sea surface temp data

Make raster with just the mean sea surface temp over the 2008-2012 time period.

Get depth data with the same resolution as the sea surface temp mean and crop to the same extent as the sea surface temp raster.

Determine which locations are suitable for Oyster aquaculture using depth and temperature parameters.

Reclassify temperature raster:

Reclassify depth raster:

Combine the reclassified depth and temp layers to find areas that are suitable for oysters:

Determine suitability in the EEZ to investigate which areas have the most potential for oyster marine aquaculture.

Suitability area visualization!

Map for the Percent Suitable Area in the Economic zones:

Map for the Total Suitable Area in the Economic zones:

Develop a function to automate the processing of determining suitable area for marine aquaculture of other species on the west coast.

The function `species_maps` :

Testing the function on a new range of temps and depths!

Oyster Aquaculture Suitability in EEZ on the West Coast

Luna Herschenfeld-Catalan

2022-12-3

Background

Goals

Datasets

1. Sea Surface Temperature

2. Bathymetry

3. Exclusive Economic Zones

Load libraries

Read in the Data

Check the crs for the data

Reproject

Prepare the depth and sea surface temp data

Make raster with just the mean sea surface temp over the 2008-2012 time period.

Get depth data with the same resolution as the sea surface temp mean and crop to the same extent as the sea surface temp raster.

Determine which locations are suitable for Oyster aquaculture using depth and temperature parameters.

Reclassify temperature raster:

Reclassify depth raster:

Combine the reclassified depth and temp layers to find areas that are suitable for oysters:

Determine suitability in the EEZ to investigate which areas have the most potential for oyster marine aquaculture.

Suitability area visualization!

Map for the Percent Suitable Area in the Economic zones:

Map for the Total Suitable Area in the Economic zones:

Develop a function to automate the processing of determining suitable area for marine aquaculture of other species on the west coast.

The function species_maps :

Testing the function on a new range of temps and depths!

The function `species_maps` :